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Abstract 

Historic events such as the uplift of mountains and climatic oscillations in the Quaternary periods greatly affected the 
evolution and modern distribution of the flora. We sequenced the tmL-tmF, ndhJ-tmL and ITS from populations throughout 
the known distributions of O. longilobus and O. taihangensis to understand the evolutionary history and the divergence 
related to the past shifts of habitats in the Taihang Mountains regions. The results showed high genetic diversity and 
pronounced genetic differentiation among the populations of the two species with a significant phylogeographical pattern 
(A/ ST >G ST , P<0.05), which imply restricted gene flow among the populations and significant geographical or environmental 
isolation. Ten chloroplast DNA (cpDNA) and eighteen nucleus ribosome DNA (nrDNA) haplotypes were identified and 
clustered into two lineages. Two corresponding refuge areas were revealed across the entire distribution ranges of O. 
longilobus and at least three refuge areas for O. taihangensis. O. longilobus underwent an evolutionary historical process of 
long-distance dispersal and colonization, whereas O. taihangensis underwent a population expansion before the main uplift 
of Taihang Mountains. The differentiation time between O. longilobus and O. taihangensis is estimated to have occurred at 
the early Pleistocene. Physiographic complexity and paleovegetation transition of Taihang Mountains mainly shaped the 
specific formation and effected the present distribution of these two species. The results therefore support the inference 
that Quaternary refugial isolation promoted allopatric speciation in Taihang Mountains. This may help to explain the 
existence of high diversity and endemism of plant species in central/northern China. 
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Introduction 

During the last three million years, repeated glacial and 
interglacial cycles have gready affected the landform and fauna 
of the total earth [1—3]. Historical ecological and biogeographic 
factors are often considered to have played significant roles in 
shaping global biodiversity by influencing regional differences in 
speciation, extinction, and migration [4-6] . The genetic structures 
of the current species record the simultaneous consequences of two 
fundamental processes: population dynamics in response to past 
geological or climatic changes, and lineage sorting within a species 
under natural selection. Thus, knowledge of the evolutionary 
history of many plant species is central to the identification of 
divergence and speciation processes [3,7-8]. 

Phylogeographical analyses can provide an understanding of 
how paleo-environmental changes in landscape and climate have 
influenced species distributions and population demography [9]. 
Such analyses play an important role in understanding the 
evolutionary history of species with changes [3]. Over the last 
decade, molecular research on the evolutionary history of plant 
species, both in China and throughout East Asia, has been 
conducted with reference to past climatic oscillations [10-18]. One 
of the results of such research is the proposal by some authors that 



glacial refugia were maintained in both the northern and the 
southern regions, or at different spatial-temporal scales in China, 
during these glacial periods. These refugia are thus suggested to 
have acted as sites for subsequent range expansion during the 
interglacial (or postglacial) periods [12-14,19-20]. Some plant 
species have extended their distribution from southwestern China 
to central/northern China [21]. Reduction of genetic diversity is 
expected under a scenario of rapid postglacial expansion, as has 
been found in northern Europe and America [2] . Because central/ 
northern China is believed not to have been glaciated during the 
Quaternary [2,22-23], it is an open question whether species 
experienced past northern expansion and southern retreat during 
the Quaternary climatic oscillations. 

Taihang Mountains locates on 35 o 19'-40°51'N and 113°10'- 
1 15°48'E. As a natural boundary mountain of the east, southeast 
of Shanxi province and Hebei and Henan provinces, it is the 
important mountain range and geographical boundary for eastern 
China. The northern part of Taihang Mountains is higher than its 
southern part. Most areas of Taihang Mountains are higher than 
1 200 meters. The uplift of Taihang Mountains, followed by the 
formation of high mountains and deep valleys within the plateau 
[24-25], was one of the most important geological events. Several 
lines of evidence suggest that the rapid uplift of Taihang 
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Mountains took place after the late Pleistocene [26-27]. Rich 
sources of species, most of which are endemic, are found in this 
region [24]. The high species richness in Taihang Mountains has 
led to the hypothesis that this region is a distribution and diversity 
center for many plant genera [24]. 

Molecular techniques have provided powerful tools for studying 
the phylogeography or migratory footprints of species [7-9,28- 
29]. Maternally inherited chloroplast DNA lineages in natural 
populations often display geographical structure [9], which is 
useful for deciphering the evolutionary history of species. The 
cpDNA markers are thought to be the most appropriate 
candidates because of their slow evolution and lack of recombi- 
nation [30]. Nevertheless, the joint use of molecular markers 
derived from different genomes provides a more complete 
description of population structure and insights into population 
history and dynamics, particularly for comparisons of maternally 
inherited organelle and bi-parentally inherited nuclear markers 
[31]. Several studies about inter-specific and intra-specific 
divergence from China, using different molecular markers, have 
played an important role in discovering phylogeographical 
patterns in East Asian flora. This is particularly true for the 
mechanisms of plant speciation, along with the production of the 
high plant biodiversity and endemism found in East Asian flora 
[12-13,15,20,32]. However, most of the studied species from 
China have been trees [13,33-34]. Herbaceous vegetation has 
received much less attention [16]. Therefore, it would be of great 
interest to study the phylogeography of an herbaceous species to 
understand the evolution and modern distribution of the 
vegetation, especially in Taihang Mountains of central/ northern 
China. 

Opisthopappus longilobus Shih and Opisthopappus taihangensis 
(Ling) Shih belong to Opisthopappus (Asteraceae). The genus 
Opisthopappus is endemic to China, and its wild distribution is 
mainly restricted to Taihang Mountains across the provinces of 
Shanxi, Hebei, and Henan, occurring on the slopes at an elevation 
of about 1000 m or in the cracks of the steep cliffs [35]. O. 
longilobus is mainly distributed in the province of Hebei, whereas 
O. taihangensis is found in the provinces of Shanxi and Henan 
[35-36]. Both species are diploid (2n= 18) [25]. A strict 
morphological differentiation occurred between O. longilobus 
and O. taihangensis. The leaf blade of O. longilobus is smooth and 
subcylindrical, and most stem and leaves are pinnatifid, with one 
pair of bracteal leaf. O. taihangensis, apprised pubescent on both 
surfaces of leaf blade, stem, and leaves, is bipinnatifid, with no 
bracteal leaf [35-36]. O. longilobus and O. taihangensis have been 
listed among the Class II State-Protected Endangered Plant 
Species [36] due to the narrow distribution range, unique 
ecological environment, and serious artificial plucking. 

O. longilobus and O. taihangensis are an ideal candidate for 
investigating the influences of Quaternary climate change in 
Taihang Mountains for the following reasons: (i) Geographically, 
O. longilobus and O. taihangensis are distributed only in Taihang 
Mountains, including Shanxi, Hebei, and Henan provinces. This 
geographic range may be advantageous for enabling investigation 
into whether the glacial refuges were maintained in Taihang 
Mountains during the last glacial maximum, or earlier cold 
periods, (ii) 0. longilobus and O. taihangensis have intra-specific 
and endemic taxa in China that display visible phenotypic 
differences, which are important in exploring the incipient 
speciation caused by isolation provided by a glacial refuge. Finally, 
(iii) 0. longilobus and 0. taihangensis are an herbaceous plant and 
sensitive to environmental changes. In research on the effects of 
climate oscillations on plant evolution, it can represent the 
phylogeography of an herb, and can serve as a model to reveal 



the historical and evolutionary processes of such plants since the 
period of Quaternary climate change. 

Thus, our objectives included the following: (1) explore the 
phylogeographical structure and the evolutionary history of O. 
longilobus and O. taihangensis; (2) elucidate the effects of the uplift 
of Taihang Mountains and the Quaternary climatic oscillations on 
the genetic structure, divergence and distribution of O. longilobus 
and O. taihangensis; (3) infer the possible refugia of 0. longilobus 
and 0. taihangensis across their whole distribution range. 

Materials and Methods 

Ethics statement 

This study was conducted in accordance with all People's 
Republic of China laws. No specific permits were required for the 
described field studies because all researchers collecting the 
samples had introduction letters from the College of Life Science, 
Shanxi Normal University, Linfen. 

Plant sampling 

Thirteen populations of Opisthopappus were investigated in this 
study, covering the entire distribution area on Taihang Mountains. 
The samples comprised five populations of O. longilobus and eight 
populations of O. taihangensis (Table 1, Figure 1). Each popula- 
tion included 10 to 25 individuals that were collected at least 5 m 
apart. Healthy leaves were collected and dried in silica gel until 
total DNA was extracted. In each of the 13 populations, the 
ecological factors of each location were determined and recorded 
(Table 1). 

DNA extraction, amplification and sequencing 

Total DNA was extracted from the silica gel-dried leaves using 
the modified 2xCTAB procedure [37]. DNA quality was checked 
by electrophoresis on 0.8% agarose gels, and DNA concentration 
was determined using an Eppendorf biophotometer protein 
nucleotide analyzer (Eppendorf China Ltd., Beijing, Germany). 
The DNA samples were diluted to 10 ng-uL -1 and stored at 20°C 
for subsequent use. 

The trnh-trnF sequences were amplified using trnh and trnF 
[38]. The primers ndhj and trnh [39] were used to amplify the 
ndhj-trnh sequences. ITS sequences were then amplified using the 
primers ITS4 and ITS5 [40]. A polymerase chain reaction (PCR) 
was then performed in a 25 |iL volume, with 50 ng plant DNA, 
2 xMasterMix (0.2 mM dNTPs, 3 mM MgCl 2 , 1 xPCR buffer, 
and 0. 1 unit Taq DNA polymerase), and 0.6 mM of both forward 
and reverse primers. The PCR parameters for all amplification 
programs of ITS and cpDNA were as follows: 4 min of pre- 
denaturation at 94°C, followed by 34 cycles of 30 s of 
denaturation, 40 s of annealing (49.2°C for ITS or 58.8°C for 
trnh-trnV or 59.8°C ndhJ-lrnL), 1 min 20 s of elongation at 72°C, 
and a final elongation step of 7 min at 72°C. PCR products were 
purified using the Wizard PCR Preps DNA purification system 
(Promega, Madison, WI, USA) following the manufacturer's 
instructions. Cycle sequencing reactions were conducted using 
the purified PCR product, AmpliTaq DNA polymerase, and 
fluorescent BigDye terminators. The sequencing products were 
analyzed using an ABI Prism 310 DNA sequencer (Applied 
Biosystems Inc., Foster City, CA, USA). 

Data analysis 

The sequences were edited manually based on the chromato- 
grams and aligned by CLUSTAL X [41] and then adjusted 
manually. Inserts and indels within all cpDNA and nrDNA 
sequences were firstiy treated as a single character resulting from 
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Figure 1. Map of sample sites for Opisthopappus (Asteraceae) on the Taihang Mountains. Location details are given in Table 1, and locality 
numbers correspond to those in Table 1 . Yellow circle dots represent the populations of 0. longilobus; red circle dots represent the populations of 0. 
taihangensis. 

doi:1 0.1 371 /journal.pone.01 04773.g001 
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one mutation. Haplotype diversity (h) and nucleotide diversity (n) 
were calculated for each population (h, n) and at the species level 
(h A , n d ) using DNAsp [42]. The program PERMUT [44] was used 
to calculate the within-population diversity (h$), total diversity (h-r), 
geographical total haplotype diversity (Ft), geographical average 
haplotype diversity (V$), level of population differentiation at the 
species level (Gst), an d an estimate of population subdivisions for 
phylogenetically ordered alleles (N ST ). We further tested the 
phylogeographical structure at the species range between (Nst) 
and (Gst) by using the U-statistic. An iVsT value higher than the 
Gst value indicates that closely related haplotypes are observed 
more often in a given geographical area than would be expected 
by chance [43]. 

To quantify the genetic differentiation partitioned among 
different groups and total genetic variance, analyses of molecular 
variance (AMOVA) were carried out using ARLEQUIN [44], 
with 1000 random permutations to test for significance of 
partitions. The spatial genetic structure of haplotypes was analyzed 
by spatial analysis of molecular variance using SAMOVA [45]. 
This program uses a simulated annealing approach to define 
groups of populations (IC) that are geographically homogenous and 
maximally differentiated from each other. In this analysis of 
haplotypes, K varied from 2 to 13, with each simulation starting 
from 100 random initial conditions. An Fqi index of genetic 
differentiation among initial K groups was computed, followed by 
an iterative simulated annealing process to obtain the optimal 
configuration of groups and final F C t values. The simulated 
annealing process was repeated 1000 times. The configuration 
with the largest F C t value among the 100 tested was retained as 
the best grouping of populations. 

Genealogical relationships between haplotypes were inferred 
from a maximum parsimony median-joining network calculated in 
NETWORK 4.5.0.2 [46]. To complement the results of 
NETWORK, TCS 1.21 [47] was also used to construct haplotype 
relationships. Phylogenetic analyses were conducted for the 
nrDNA and cpDNA sequence data using maximum likelihood 
(ML) and Bayesian inference (BI), respectively. 

The time of the most recent common ancestor (TMRC A) of all 
haplotypes was estimated via a Bayesian approach implemented in 
*BEAST (Star BEAST) [48-49] using a GTR+I substitution 
model. The best substitution model was determined according to 
the Akaike Information Criterion (AIC) in jModeltest [50] . These 
models were applied in ML, BI, and subsequent *BEAST analyses. 
Bayesian inference was performed using MrBayes 3.1.2 [51]. Two 
independent runs of Metropolis-coupled Markov chain (MCMC) 
analysis were executed, each including one cold chain and three 
incrementally heated chains that started randomly in the 
parameter space. Two independent runs of 10 generations were 
carried out, with sampling at every 1000 generations. The first 
25% of sampled trees were discarded as burn-in, and the 
remaining trees were used to construct a Bayesian consensus tree. 
The convergence of chains was checked using Tracer 1.5 [52]. 
The remaining trees were pooled to estimate the posterior 
probabilities (PPs). 

No fossil records or biogeographic events isolating distinct 
populations are available to calibrate a cpDNA-IGS substitution 
rate for 0. longilobus and O. taihangensis. Therefore, for our 
combined chloroplast non-coding regions, we assumed minimum 
and maximum values of a range of average mutation rates 
reported for synonymous sites of plant chloroplast genes [i.e., 1 .2 
and 1.7x10 9 substitutions per site per year (s/s/y)] [53]. These 
rates were then used for estimating TMCRA in *BEAST under a 
relaxed molecular clock assumption. 



To detect whether population groups experienced recent 
population expansion and satisfied the assumption of neutrality, 
a mismatch distribution analysis was performed using the DNAsp 
program. Tajima's D and Fu and Li's D* were calculated for the 
entire genus and groups of populations. The statistical significance 
of D and D* was estimated with coalescent simulations as 
implemented in this program [54—55]. To further infer demo- 
graphic processes, the null hypotheses of a spatial expansion and of 
a pure demographic expansion were tested in ARLEQUIN by 
comparing observed and expected distributions of pairwise 
sequence differences (mismatch distributions). For each expansion 
model, goodness-of-fit was tested using the sum of squared 
differences (SSD) and Harpending's [56] raggedness index 
(HRag). 

Finally, isolation by distance (IBD) for both nrDNA and cpDNA 
data was tested based on pairwise geographical and genetic 
distances (F ST ) with the Isolation by Distance Web Service [57] 
(http://ibdws.sdsu.edu/~ibdws/), running 10,000 Mantel per- 
mutations. 

Results 

Patterns of variability in cpDNA and ITS sequences 

The trnh-lrnF and ndhj-trnh intergenic spacers varied from 
800 to 866 bp and from 790 to 835 bp, respectively. The trnh- 
trn¥ intergenic spacer, which exhibits considerable nucleotide 
polymorphism, is characterized by 10 substitutions and a 4 bp 
insertions and deletions (Table SI). The 4 bp insertions and 
deletions region was mainly observed in 0. longilobus populations. 
The ndhJ-trnL. spacer was characterized by one indel, three 
substitutions, and a 10 bp insertions and deletions (Table SI). The 
total cpDNA combined matrix comprised 1701 sites, of which 13 
positions were variable and 18 harbored gaps. Ten haplotypes 
(Figure 2A) were observed by combining cpDNA data: eight 
haplotypes were within O. longilobus populations and five were 
within 0. taihangensis populations. The haplotype diversity (h) 
ranged from 0 to 0.854, and the nucleotide diversity (n) from 0 to 
0.00194 in the 5 populations of O. longilobus (Table 1). At species 
level, h d = 0.838 and n d = 0.00 175 for 0. longilobus. Within 8 0. 
taihangensis populations, h ranged from 0 to 0.667, n from 0 to 
0.00107, h d was 0.562 and n d was 0.00073. The highest haplotype 
diversity was observed in population WDS from 0. longilobus and 
BJY from O. taihangensis (Table 1). 

The full ITS sequences, including ITS1+5.8S+ITS2, were 
682 bp in length, comprising of 254 bp and 224 bp for ITS1 and 
ITS2, respectively. Eighteen haplotypes (Figure 2B) and 17 
polymorphic sites, the latter consisting of 15 parsimony informa- 
tive sites and 2 singleton variable sites, were detected at ITS 
sequences. Two populations, SBY and XTS, were monomorphic, 
whereas the remaining populations were polymorphic (Table 1). 
The haplotype diversity (h) ranged from 0 to 0.857 and the 
nucleotide diversity (n) from 0 to 0.00391 in the 5 populations 0. 
longilobus. Within 8 0. taihangensis populations, the haplotype 
diversity h ranged from 0 to 0.714 and the nucleotide diversity n 
from 0 to 0.00293. For 0. longilobus, h A = 0.908, n d = 0.00537; for 
O. taihangensis populations, h d = 0.558, n d = 0.00158 (Table 1). 

Genetic diversity and population structure 

A genetic diversity analysis revealed high genetic diversity 
harbored in O. longilobus and 0. taihangensis populations. The 
parameters h T , Ft, As, and V$ were 0.908, 0.916, 0.669, and 
0.581, respectively, for O. longilobus based on cpDNA sequences. 
By ITS sequences, At = 0.944, Ft = 0.937, h s = 0.824, and 
V s = 0.889 within 0. longilobus populations. For O. taihangensis 
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Figure 2. Geographical distribution of 10 cpDNA (A) and 18 (B) nrDNA haplotypes. The pie charts reflect the frequency of haplotype 
occurrence in each population. Haplotype colours correspond to those shown in panel. 
doi:1 0.1 371 /journal.pone.01 04773.g002 



populations, h T = 0.828, V T = 0.833, h s = 0.531, and F s = 0.455 
by cpDNA sequences; h T = 0.806, Ft — 0.810, h s = 0.704, and 
Fs = 0.592 with nrDNA sequences. 

Based on ITS sequences, phylogenetic analyses were carried 
(Figure SI). All studied populations made up a single monophyletic 
lineage, indicating that all of these divergent populations were 
derived within Opisthopappus itself, rather than the result of 
introgression by other species. Spatial genetic analyses of cpDNA 
and nrDNA haplotypes in all 13 populations using SAMOVA 
indicated that Fq-y increased to a maximal value of 0.921 when 
K — 2 (K, the number of groups). Thus, the division of all the 13 
sampled populations approximately into two groups is appropri- 
ate. 

AMOVA analysis revealed that 60.58% and 67.70% (P<0.01) 
of the total variation was due to differences among the cpDNA 
and nrDNA populations for all studied populations, respectively 
(Table 2). When the populations were grouped by taxonomic 
variety, the cpDNA showed 27.28% variation among populations 
within species (P<0.01; Table 2). Similarly, for the nrDNA, 
AMOVA revealed that 26.68% of the variation was attributed to 
differences among populations within the species (P<0.01; 
Table 2). There was a significant (P<0.01) variation between 
species, suggesting that the two separately distributed groups have 
a genetic differentiation. In addition, 60.71% of the total cpDNA 
variation and 68.24% of the total nrDNA variation existed among 
O. longilobus populations (Table 2). The genetic structure of 0. 
taihangensis populations showed a similar trend with that of 0. 
longilobus populations, with most of the genetic differentiation 
occurring among populations (Table 2). 



Results of the Mantel test showed a significant correlation 
between nrDNA data differentiation among populations and the 
natural logarithm of the geographical distances throughout the 
sampled range of O. longilobus (r = 0.442; P = 0.001) and O. 
taihangensis (r = 0.385; P = 0.001). Likewise, a significant corre- 
lation was detected among populations from cpDNA data of O. 
longilobus (r = 0.401; P — 0.001) and 0. taihangensis (r = 0.399; 
P = 0.001). This indicates a significant correlation between genetic 
differentiation and geographical distance, implying strong IBD. 

Haplotype and phylogeographical structure 

In this study, similar topologies were obtained from the two 
different network approaches used to infer the relationships among 
the cpDNA and nrDNA haplotypes. Only the median-joining 
network is shown. The haplotype H8 was the most frequent 
haplotype in cpDNA data, occurring in 6 of the 13 populations. 
Three haplotypes, HI, H3, and H8, were shared by 0. longilobus 
and 0. taihangensis populations (Table 1 and Figure 3). In 
nrDNA data, Hnl was the most frequent haplotype, occurring 
in 10 of the 13 populations (Figure 3 and Figure 4). The 
haplotypes Hnl, Hn4, and Hn5 were shared by 0. longilobus 
and 0. taihangensis populations (Table 1 and Figure 4). The 
distributions of haplotypes that were restricted to a single 
population were extremely skewed. Four cpDNA and nine nrDNA 
haplotypes were found in O. longilobus single populations. Only 
one cpDNA haplotype was identified for 0. taihangensis 
populations. 

Based on cpDNA data, we estimated a G ST value (0.649, P< 
0.05), which is significandy smaller than the Nst value (0.761). 
When the nrDNA was examined, the Gst value (0.253) also 
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Table 2. Analyses of Molecular Variance (AMOVA) based on the nrDNA (ITS) and cpDNA sequences. 





Source of variation 


d.f 


Sum of squares 


Variance components 


Variation percentage 


Fixation indices 


ITS{all population) 


Among populations 


19 


105.101 


1.218 


67.70 


F ST = 0.677 (P<0.01) 


Within populations 


178 


44.179 


0.581 


32.30 




ITS (species) 


Among species 


1 


54.511 


1.113 


48.16 


F CT = 0.482(P<0.01) 


Among populations within species 


18 


50.590 


0.617 


26.68 


F ST = 0.748(P<0.01) 


Within populations 


178 


44.179 


0.581 


25.16 


F SC = 0.515(P<0.01) 


ITS(0. longilobus) 


Among populations 


12 


44.767 


1.499 


68.24 


F ST = 0.682 (P<0.01) 


Within populations 


76 


23.022 


0.698 


31.76 




ITS(0. taihangensis) 


Among populations 


14 


21.157 


0.492 


90.04 


F ST = 0.900 (P<0.01) 


Within populations 


95 


5.823 


0.054 


9.96 




ITS (SMOVA species) 


Among species 


7 


57.023 


2.456 


50.21 


F CT = 0.502(P<0.01) 


Among populations within species 


12 


55.147 


2.078 


28.34 


F ST = 0.786(P<0.01) 


Within populations 


178 


40.563 


1.789 


21.45 


F SC = 0.497(P<0.01) 


cpDNA(all population) 


Among populations 


19 


78.452 


2.187 


60.58 


F ST = 0.606 (P<0.01) 


Within populations 


178 


50.234 


1.023 


39.42 




cpDNA (species) 


Among species 


1 


60.231 


1.586 


46.59 


F CT = 0.466(P<0.01) 


Among populations within species 


18 


58.427 


0.789 


27.28 


F ST = 0.738(P<0.01) 


Within populations 


178 


45.126 


0.567 


26.13 


F SC = 0.533(P<0.01) 


cpDNA (0. longilobus) 


Among populations 


12 


29.313 


1.032 


60.71 


F ST = 0.607 (P<0.01) 


Within populations 


76 


19.364 


0.668 


39.29 




cpDNA (0. taihangensis) 


Among populations 


14 


21.601 


0.467 


68.64 


F ST = 0.686 (P<0.01) 


Within populations 


95 


9.183 


0.214 


31.36 




cpDNA (SMOVA species) 


Among species 


7 


56.273 


1.756 


52.34 


F CT = 0.523(P<0.01) 


Among populations within species 


12 


55.012 


1.111 


27.55 


F ST = 0.799(P<0.01) 


Within populations 


178 


41 .603 


0.654 


22.36 


F SC = 0.498(P<0.01) 



doi:1 0.1 371 /journal.pone.01 04773.t002 



- N^ipunanthemum 



-Euppofytia 



H2Q 



o 



H10O H8 



H4 0" 

H6 0 



H5 0 



H9(0 



Figure 3. The evolutionary relationships among cpDNA haplotypes. (A) NJ phylogenetic tree for the 10 cpDNA haplotypes. Numbers above 
branches are support values from bootstrap resampling/Bayesian inference. (B) Median-joining network. Sizes of the circles are proportional to the 
overall frequency of the haplotypes in the entire sample of all species. 
doi:1 0.1 371 /journal.pone.01 04773.g003 
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Figure 4. The evolutionary relationships among nrDNA haplotypes. (A) NJ phylogenetic tree for the 18 nrDNA haplotypes. Numbers above 
branches are support values from bootstrap resampling. (B) Median-joining network. Sizes of the circles are proportional to the overall frequency of 
the haplotypes in the entire sample of all species. 
doi:1 0.1 371 /journal. pone.0104773.g004 



differed significantly from the ATst value (0.641). A significant 
phylogeographic structure was indicated by both cpDNA and 
nrDNA data (JV ST >G ST , P<0.05). 

The phylogeographical relationships of the 1 0 cpDNA haplo- 
types and 18 nrDNA haplotypes were assessed under Neighbor- 
joining analyses (NJ) and Bayesian inferences drawn using 
Hippolytia and Nipponanthemum as outgroups (Figure 3 and 
Figure 4). Similar topology to phylogeographical relationships 
indicated a single phylogenetic split. All haplotypes were clustered 
into two lineages: I and II (Figure 3 and Figure 4). Lineage I 
contained some populations of 0. longilobus with a high bootstrap 
support. The remaining haplotypes formed lineage II, containing 
some 0. longilobus populations and all 0. taihangensis popula- 
tions, also with a high bootstrap support. This unrooted network of 
cpDNA haplotypes is consistent with the strict consensus tree 
produced by NJ and Bayesian inference and furthermore displays 
the relationship of interior (ancestral) and tip (derived) haplotypes 
(Figure 3). For lineage I (Figure 3), the haplotype H3 located in 
the center and fixed in three populations, being an ancestral 
haplotype. In clade II, HI was shared by five populations. It 
appears to be an ancestral haplotype distributed in the center of 
clade II. The haplotype H8 was located in the tip of the network 
despite though contained a higher frequency of haplotypes. In the 
nrDNA haplotype network (Figure 4), Hn8 was ancestral for the 
lineage I; Hnl was the predominant and widespread one in the 
clade II. 

To test whether the current expansion had occurred in different 
populations, we calculated the frequency distribution of pairwise 
nucleotide differences among 0. longilobus populations and 0. 
taihangensis populations respectively. The resultant mismatch 
distribution consisted of multiple multimodal curves (Figure S2). 
Neutrality tests revealed nonsignificant positive values in consid- 
ering all populations (Table 3). The observed mismatch distribu- 
tions of cpDNA haplotypes were multimodal for all groups except 
for O. taihangensis. The mismatch distributions based on nrDNA 
haplotypes were present a similar pattern with cpDNA data 
(Table 3). Thus, a recent population expansion for 0. longilobus is 
unlikely. 0. taihangensis perfectly fits the expected expansion 
distribution (Figure 3). 

Divergence times were estimated using the program *BEAST. 
The ESS values were from 150 to 240 for all the nodes discussed 
as follows. All sampled haplotypes of Opisthopappus coalesced at 



about 1.90 Ma (95% highest posterior density, HPD) (0.50-1.72), 
which indicated the origin of Opisthopappus during the late 
Tertiary periods. The divergence times, estimated at approxi- 
mately 1.28 Ma (95% HPD) (0.70-2.34) between 0. longilobus 
and 0. taihangensis , are shown in the haplotype phylogenic tree 
(Figure 3A). O. longilobus and 0. taihangensis thus appear to have 
diverged from each other during the early Pleistocene. 

Discussion 

Taxonomic implications 

At the morphological character levels, the evidences that 
illustrated O. taihangensis as a separate species from O. longilobus 
were only leaf form and involucres. So, these two species were 
once regarded as a species. In this study, the analyses of the 
combined cpDNA and nrDNA matrix support the monophyly of 
two clades (Figure 3A, Figure 4A, and Figure SI). The phyloge- 
netic relationships among all populations show an overall 
congruence with the result of the haplotype network (Figure 3B 
and Figure 4B). All haplotypes are split into two major groups 
according to cpDNA and nrDNA data. The populations from two 
distinct clades approximately exhibited specific geographical 
distributions in Taihang Mountains. SAMOVA analysis identified 
two approximately defined groups corresponding to O. longilobus 
and O. taihangensis populations at every level of divergence. 
Within DNA sequences of cpDNA and ITS, the total of thirty 
polymorphic sites harbored in populations. Unique haplotypes 
found in O. longilobus and O. taihangensis populations respec- 
tively. Moreover, a significant variation based on AMOVA 
occurred between 0. longilobus and 0. taihangensis 
(^stc P dna = 0.482 and F STnrDNA = 0.466, Table 2). By used 
SPvAP and cpSSR markers [58-59], AMOVA indicated that a 
significant genetic differentiation between O. longilobus and 0. 
taihangensis. Accordingly, our molecular results support that 0. 
longilobus and 0. taihangensis were regarded as two species rather 
than a species. That result is consistent with the point of Hu [25]. 
The divergence times between O. longilobus and 0. taihangensis 
(Figure 3A) indicate that genetic divergence between two species 
occurred in early Pleistocene. This split was related to the 
ecological and geographical habitat changes resulting from climate 
oscillations during the Quaternary glacial periods and complicated 
topography with the uplift of Taihang Mountains. 
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Geographic patterns of genetic diversity and structure 

The high levels of genetic diversity were observed in O. 
longilobus and O. taihangensis populations (Table 1). These 
results are similar to previous research on herbaceous plants in 
East Asia and significandy exceed the average value of 0.67 for the 
170 plant species documented [8,12,19,31,60], In O. longilobus, 
the WDS population obtained the most haplotypes and highest 
genetic diversity, followed by the SHS population. This suggests 
the presence of two genetic diversity centers for O. longilobus, i.e., 
the regions where WDS and SHS populations are located. For O. 
taihangensis, the BJY and WML populations, which obtained the 
high genetic diversity, were genetic diversity centers for O. 
taihangensis. Significant genetic divergence and a highly structure 
were illustrated in all studied populations (Figure 3, Figure 4, and 
Table 2). For O. longilobus (i^sTcpDNA = 0.607, -FsTm-DNA = 0.682) 
and O. taihangensis (^sTcpDNA - 0.686, F STnrWA = 0.805), the 
results of our analysis suggested that significant variation was 
detected from among populations rather than from variation 
within groups (Table 2). Moreover, very high genetic differenti- 
ation in O. longilobus and O. taihangensis greatly exceeding the 
mean value (G ST = 0.22) reported by Nybom [61]. 

Species genetic diversity and structure can be affected by life 
history traits (e.g., life cycle, breeding system, pollination 
mechanism) and environmental effects (e.g., geographical range, 
climate, topography) [61-62]. O. longilobus and O. taihangensis 
possess a sexual reproductive mode and a relatively long life cycle 
as perennial herbs, all of which are traits associated with low total 
genetic diversity [63-64]. However, high genetic diversity was 
observed in our study, and this diversity appears to result from 
high among-population genetic variation. It has been widely 
documented that the dispersal of pollen and seeds in plant species 
is strongly linked to the development of a population genetic 
structure in bi-parentally inherited (nuclear-encoded) and mater- 
nally inherited (cytoplasmic) genetic markers [12,31,61,62,64]. 
Generally, outcrossing taxa with high seed dispersal capacity 
retain the majority of both types of genetic marker variation within 
populations, whereas selfing (and/ or asexual) taxa with restricted 
seed dispersal allocate the majority of such variation among 
populations. O. longilobus and O. taihangensis are self-compatible 
plant species [25] . These two species, therefore, might be expected 
to have a heterogeneous distribution of both cpDNA and nrDNA 
variation among populations. Taken at face value, the above 
results may suggest that restricted gene flow via pollen and seed 
has resulted in significant population genetic differentiation in O. 
longilobus (A r ml . pD NA = 0.12; A r mnrDNA = 0.10) and O. taihangensis 
(NmcpDNA- 0.19; W mnrD NA = 0.14). The physiographic complexity 
of Taihang Mountains and the deeply carved valleys and ravines 
between the inhabit areas probably would imposes significant 
barriers to gene flow among populations of O. longilobus and O. 
taihangensis. The genetic divergence between the two lineages 
occurred at 1.28 Ma, corresponding to the early Pleistocene. This 
divergence falls within the glaciations during the Taku Glaciation 
(1.05-1.20 Ma). The glaciers and/or extremely low temperature 
in the mountains might also have created barriers to gene flow 
between geographically isolated populations, which therefore 
promoted the divergence. For both marker systems, we generally 
observed a good fit of these populations to an IBD model. When 
coupled with high levels of population subdivision, such conditions 
strongly suggest that geographical isolation had a larger historical 
role in extant population structure compared with limited pollen 
and seed flow alone. The habitats of O. longilobus and O. 
taihangensis are either rocky or massifs, and these species grow 
discontinuously in different habitats along the altitudinal gradient 
from 700 to 1500 m, in which the variable micro-surrounding, 
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complex topography, and great altitudinal variability of the region 
might promote the high degree of genetic population differenti- 
ation [65-66]. Limited migration coupled with a heterogeneous 
environment could promote local adaptation and fixation of 
different alleles in O. longilobus and O. taihangensis populations. 

Phylogeographical structure and inference of 
demographic history 

Both 0. longilobus and 0. taihangensis populations have a 
significant phylogeographic structure. The current distribution of 
haplotypes could be a result of past climatic changes related to the 
advance/retreat of glaciations [11]. During the Quaternary 
periods, climatic oscillations resulted in repeated drastic environ- 
mental changes, which further caused massive range shifts of most 
plants and animals, leading to accumulated genetic differences and 
particular phylogeographic patterns [2-3] . 

Our analyses suggested that Opisthopappus originated during 
the Pliocene. Before the Pleistocene, the current range of Taihang 
Mountains would have been covered by grassland with roughly 
similar geographic conditions [67]. As herbs, 0. longilobus and 0. 
taihangensis can disperse following grassland migration. On this 
basis, we suspect that 0. longilobus and O. taihangensis were once 
distributed throughout Taihang Mountains and its adjacent 
region. With uplift of Taihang Mountains [27] and climate 
fluctuation of Quaternary glacial periods, 0. longilobus and 0. 
taihangensis were fragmented and forced into the refuge areas. 
The flora would migrate and expand during interglacial and 
glacial periods. This migration and expansion resulted in 
ecological displacement, which might have resulted in populations 
situated at different spatial-temporal scales. Different habitats 
might enhance the isolation of plant populations and prevent gene 
flow, which in turn could lead to incipient allopatric speciation or 
further allopatric speciation [12,18,19,68,69]. During the Quater- 
nary periods, the paleovegetation in Taihang Mountains repeat- 
edly appear the replacement of grasslands and forests cycles [67]. 
Those multiple replacements likely fragmented and isolated the 
habitat of herbs. Therefore, the refuge areas for 0. longilobus and 
O. taihangensis would have remained separated and fragmented at 
different spatial-temporal scales for a long time. Moreover, 
dramatic geological changes induced by the Taihang Mountains 
uplift contributed to the extremely complicated landscape, which 
likely acted as a major driving force for the separation of the 0. 
longilobus and O. taihangensis. Divergent selection between 
populations in contrasting environments, longer temporal and 
spatial separation and fragmentation would accelerate their 
differentiation, and eventually promote allopatric speciation. 

Following this divergence, 0. longilobus and 0. taihangensis 
underwent strikingly different histories. For O. longilobus, long- 
distance dispersal and colonization may be a major historical 
process. 0. longilobus was divided into different lineages around 
0.89 Ma, which fell within the interglacial stage between the Taku 
Glaciation and Lushan Glaciation. The paleovegetation grasslands 
of Taihang Mountains [67] might be regarded as a corridor 
connecting the currently isolated populations. Populations of 0. 
longilobus in Hebei province may have extended across this 
corridor and reached Henan province. Haplotypes HI and Hnl 
located in populations LLS and WDS or LLS and SBY imply this 
long dispersal (Table 2). The 10 bp insertions of ndhj-trnh 
occurring in population LLS also supported this viewpoint (Table 
SI). Subsequently, with the advent of inter-glaciations or post- 
glaciations, the climate changed, which led to the paleovegetation 
of Taihang Mountains experiencing a transition from grasslands to 
forests. The corridor for O. longilobus migration gradually 
disappeared, which would have resulted in isolation of the 



population [11,12,68]. Additionally, we did not detect any 
population expansion due to the clear multimodel mismatch 
distribution (Figure S2). The complex geological conditions of 
Taihang Mountains would limit 0. longilobus population expan- 
sion on a macro scale and thus restricting O. longilobus to sustain 
in situ during the Quaternary periods, perhaps by moving 
upwards and downwards in their mountain ranges. For 0. 
taihangensis, a pattern of population expansion was observed, 
because the observed mismatch distribution for cpDNA and 
nrDNA haplotypes was unimodal (a pattern consistent with 
expansion) (Figure S2). This expansion was dated to be 
0.78 Ma, which corresponding the interglacial stage between the 
Taku Glaciation and Lushan Glaciation. The above results suggest 
an expansion consistent with a range expansion and re-coloniza- 
tion before the mainly uplift of Taihang Mountains. With the 
uplift of Taihang Mountains, large geographical distances and 
significant topographical barriers could restrict gene flow and limit 
any large-scale population expansion [12,70]. The cpDNA 
sequences characteristic of the trnh-lrnF and ndh]-trnL intergenic 
spacers further confirmed this hypothesis (Table SI). 

Crandall and Templeton [7 1] suggest that ancestral haplotypes 
are located in a central position within the haplotype phylogeo- 
graphical network and that potential refugia are usually charac- 
terized by widely dispersed ancestral haplotypes and other tip and 
unique haplotypes. Geographic areas displaying increased levels of 
genetic diversity are thus good candidates in the search for past 
refuges [8,34] . These regions are characterized by relatively stable 
ecological conditions during environmental fluctuations, which 
foster the accumulation of genetic diversity [72]. Most 0. 
longilobus individuals or populations gathered in lineage I 
(Figure 3 and Figure 4). Haplotypes H3 and Hn8 are located in 
the center of the haplotype network. And Haplotype H3 is 
detected in three populations (WDS, SHS and WML), Hn8 in two 
populations (WDS and SHS) (Figure 2). WDS and SHS are 
located in the southern part of Hebei province and have higher 
haplotype diversity (Table 1). Thus, we suggest that the WDS and 
SHS populations' regions were two major potential refugia for O. 
taihangensis. These regions may therefore have been character- 
ized by a relatively mild climate, potentially containing microcli- 
matic environments able to impart relative stability to a range of 
habitats. Lineage II consisted of all populations of O. taihangensis 
and some individuals of 0. longilobus. Haplotypes HI and Hnl, 
which are located in the center of the haplotype network, were the 
most dominant haplotypes found in five and ten populations, 
respectively. Moreover, some populations in these haplotypes (HI 
and Hnl) had higher haplotype diversity (Table 1, e.g., BJY, 
WML, and YTS). Thus, we suggest that there were multiple small 
refuges located in Shanxi and Henan provinces during the last 
glacial maximum or earlier cold periods for 0. taihangensis. The 
above populations of 0. taihangensis, distributed in southern 
Taihang Mountains, would have been provided with a complex 
and stable environment. Thus, the potential refugia area may have 
been located in those regions. 

However, our comprehensive sampling showed that assessing 
the species-level relationships in our study is complicated by the 
fact that not all cpDNA and nrDNA haplotypes are species- 
specific. Incongruence between the species boundaries and the 
genealogy of their cpDNA and nrDNA sequences is illustrated by 
the sharing of haplotypes (e.g., HI and Hn8) between the two 
species. On the other hand, several haplotypes were found within a 
single species (Table 1 and Figure 2). Haplotype distribution does 
not strictly follow species circumscription. Moreover, an associa- 
tion of haplotypes with geographically circumscribed regions 
rather than with taxonomic boundaries is a phenomenon observed 
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in 0. longilobus and 0. taihangensis . The sharing of haplotypes 
and lack of reciprocal monophyly might be explained by the 
persistence of ancestral polymorphisms during speciation events 
and/ or exchange of genes by inter-specific hybridization (Table S 1 
and Figure SI). These two species have similar habitats and the 
same period of flowering [25], and possible visitation by the same 
pollinators in the absence of physiological or genetic reproductive 
barriers might have enabled hybridization events. Nonetheless, 
further work on chromosomes and the sequencing of additional 
genomic regions are needed to make more in-depth conclusions 
about the hybrid status of these taxa. 

Conclusions 

In summary, O. longilobus and 0. taihangensis have a high level 
of genetic diversity. The following phylogeographical structure 
and historical scenario of 0. longilobus and O. taihangensis can 
also be drawn. The vicariance following the uplift of Taihang 
Mountains and a transition of paleovegetation of Taihang 
Mountains mainly shaped the present distribution of these two 
species. The above findings imply that multiple small refugia could 
have been maintained in the Taihang Mountains regions, which 
allowed 0. longilobus and 0. taihangensis to persist in situ and 
maintain sizable haplotype and nucleotide diversity at the lineage- 
wide scale during the glaciations. The current distribution pattern 
of O. longilobus and O. taihangensis refuges is similar to that of 
other East Asian plants (e.g., Cathaya argyrophylla, Picea 
crassifolia, Sinopodophyllum hexandrum, Lagochilus Bunge ex 
Bentham) [8,20,73], with multiple refuges sustained across their 
distribution ranges. Furthermore, complex topography rendered 
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